Use dependent paper

rm(list = ls())
library(tidyverse)

Warning: package ‘tidyverse’ was built under R version 3.6.2
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.1
Warning: package ‘ggplot2’ was built under R version 3.6.2
Warning: package ‘tibble’ was built under R version 3.6.2
Warning: package ‘tidyr’ was built under R version 3.6.2
Warning: package ‘readr’ was built under R version 3.6.2
Warning: package ‘purrr’ was built under R version 3.6.2
Warning: package ‘dplyr’ was built under R version 3.6.2
Warning: package ‘forcats’ was built under R version 3.6.2
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library('svglite')

Warning: package ‘svglite’ was built under R version 3.6.2

library(broom)

Warning: package ‘broom’ was built under R version 3.6.2

library(purrr)
library(lemon)

Warning: package ‘lemon’ was built under R version 3.6.2

Attaching package: ‘lemon’

The following object is masked from ‘package:purrr’:

    %||%

The following objects are masked from ‘package:ggplot2’:

    CoordCartesian, element_render

library(pracma)

Attaching package: ‘pracma’

The following object is masked from ‘package:purrr’:

    cross

library(lubridate)

Warning: package ‘lubridate’ was built under R version 3.6.2

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

library(plotrix)

Warning: package ‘plotrix’ was built under R version 3.6.2

library(lme4)

Warning: package ‘lme4’ was built under R version 3.6.2
Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:pracma’:

    expm, lu, tril, triu

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(lmerTest)

Warning: package ‘lmerTest’ was built under R version 3.6.2

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:pracma’:

    rand

The following object is masked from ‘package:stats’:

    step

library(car)

Loading required package: carData
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'car':
  method                          from
  influence.merMod                lme4
  cooks.distance.influence.merMod lme4
  dfbeta.influence.merMod         lme4
  dfbetas.influence.merMod        lme4

Attaching package: ‘car’

The following object is masked from ‘package:pracma’:

    logit

The following object is masked from ‘package:dplyr’:

    recode

The following object is masked from ‘package:purrr’:

    some

library(sjstats)

Warning: package ‘sjstats’ was built under R version 3.6.2

Attaching package: ‘sjstats’

The following object is masked from ‘package:broom’:

    bootstrap

call_aesthethics <- function(text_size){
  
  th <- theme(   panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(),
                 panel.border = element_blank(),
                 panel.background = element_blank(),
                 axis.line = element_line(size = 0.5), 
                 legend.position = 'right', 
                 legend.text = element_text(size= text_size, family="Helvetica"),
                 text = element_text(size= text_size, family="Helvetica"), 
                 strip.text.x = element_text(size = rel(0.90)), 
                 strip.text.y = element_text(size = rel(0.90)), 
                 axis.title.x = element_text(vjust=-0.3), 
                 plot.title = element_text(hjust = 0.5, vjust = 0), 
                 axis.ticks = element_line(size = 0.4), 
                 axis.text.x.bottom  = element_text(size = rel(0.90), margin = unit(c(t = 2.5, r = 0, b = 0, l = 0), "mm")),
                 axis.title.y = element_text(vjust = 1),
                 axis.text.y = element_text(size = rel(0.90), margin = unit(c(t = 0, r = 2.5, b = 0, l = 0), "mm")),
                 axis.ticks.length = unit(-1.2, "mm"),
                 axis.text.x.top = element_text(size = rel(0.90), margin = unit(c(t = 0, r = 0, b = 2.5, l = 0), "mm")))
  
  return(th)
}

my.colors2 <-  c("#006600", "#800080", "#FF9900", 'deepskyblue4')
th <- call_aesthethics(16)

# baseline subtraction function
baseline_subtract <- function(my.data, base_cn_start, base_cn_end, returnAngle){
  
  
  my.subjects <- unique(my.data$SN)
  counter <- 1
  mybaseline <- setNames ( data.frame(matrix(NaN, nrow = length(my.subjects) * 12, ncol = 5)), 
                           c('SN', 'ti', 'Hand_base', 'RT_base','Return_base'))
                                                    
  for(si in my.subjects){
    
    idx <- my.data$SN == si
      
    for(mi in unique(my.data$ti[idx])){
      
      idx_sub <- my.data$SN == si & my.data$ti == mi
      idx_bl <- my.data$SN == si & my.data$ti == mi & 
        my.data$CN >= base_cn_start & my.data$CN <= base_cn_end
      
      sub_Hand_mean <- mean(my.data$Hand[idx_bl], na.rm = TRUE)
      my.data$Handb[idx_sub] <- my.data$Hand[idx_sub] - sub_Hand_mean
      
      sub_RT_mean <- mean(my.data$RT[idx_bl], na.rm = TRUE)
      my.data$RTb[idx_sub] <- my.data$RT[idx_sub] - sub_RT_mean
      
      if(returnAngle == 1){
        sub_Return_mean <- mean(my.data$ReturnAngle[idx_bl], na.rm = TRUE)
        my.data$ReturnAngleb[idx_sub] <- my.data$ReturnAngle[idx_sub] - sub_Return_mean
        mybaseline$Return_base[counter] <- sub_Return_mean
      }
      
      mybaseline$SN[counter] <- si
      mybaseline$ti[counter] <- mi
      mybaseline$Hand_base[counter] <- sub_Hand_mean
      mybaseline$RT_base[counter] <- sub_RT_mean
      counter <- counter + 1
    }
  }
  
  return(list(my.data, mybaseline))
}

mainDir <- "~/Dropbox/VICE/JT/KIM_USEDEPENDENT/Figures"
subDir <- today()
figDir <- "~/Dropbox/VICE/JT/KIM_USEDEPENDENT/Figures/"
knitr::opts_knit$set(root.dir = figDir)
dir.create(file.path(mainDir, subDir), showWarnings = FALSE)

Renalaysis of V&S

VS.data <- as_tibble( read.csv("/Users/jonathantsay/Dropbox/MOTOR/use_dependent/use-dependent-learning/VerstynenSabes2011.csv", header = TRUE, sep = ',')) %>% mutate(Bias = Hand_IB)

VS Individual data

for(si in unique(VS.data$SN)){
  
  ind.plot <- VS.data %>% 
    filter(SN == si & Distance > 0) %>%  
    ggplot(aes(x = RT, y = Bias, group = factor(Distance))) + 
    geom_point(alpha = 0.4, aes(color = factor(Distance))) + 
    geom_smooth(method = 'lm', aes(color = factor(Distance)), size = 0.5, se = FALSE) + 
    #stat_cor(method = 'pearson')+ 
    facet_rep_wrap(.~Distance, repeat.tick.labels = TRUE) + 
    scale_color_manual(values = my.colors2[c(1, 2, 3)], guide = FALSE) + 
    labs(y = 'Inward Bias (°)', x = 'RT (s)') + #subtitle = sprintf('V&S Subject #%s', si) 
    th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
               legend.position = 'none') +
    scale_x_continuous(breaks = c(0.2, 0.4)) + 
    scale_y_continuous(breaks = seq(-50, 100, 50)) + 
    coord_capped_cart(ylim = c(-50, 100)) 
  print(ind.plot) 
  
  ggsave(sprintf('VS_SUB%s_notitle_%s.pdf', si,subDir), plot = ind.plot, height = 3, width = 6,  units = "in")
  
}

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (geom_point).
`geom_smooth()` using formula 'y ~ x'
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (geom_point).
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

VS.group data

VS.data.ind <- VS.data %>%
  group_by(SN, Distance) %>%
  dplyr::summarise(bias_mean = mean(Bias, na.rm = TRUE), 
                   rt_mean = mean(RT, na.rm = TRUE))

`summarise()` has grouped output by 'SN'. You can override using the `.groups` argument.

VS.data.grp <- VS.data.ind %>%
  group_by(Distance) %>%
  dplyr::summarise(bias_grp_mean = mean(bias_mean, na.rm = TRUE), 
                   bias_stderr_mean = std.error(bias_mean, na.rm = TRUE), 
                   rt_grp_mean = mean(rt_mean, na.rm = TRUE), 
                   rt_stderr_mean = std.error(rt_mean, na.rm = TRUE)) %>%
  add_row(Distance = 0,  bias_grp_mean= 0, bias_stderr_mean = 0)

VS.grp <- VS.data.grp %>%
  mutate(groupvar = 1) %>%
  ggplot(aes(x = Distance, y = bias_grp_mean, group = groupvar)) + 
  geom_segment(aes(x=-5, xend=95, y = 0, yend = 0), color = 'lightgrey') +
  geom_point() + 
  geom_line() + 
  geom_errorbar(aes(ymin = bias_grp_mean - bias_stderr_mean, ymax = bias_grp_mean + bias_stderr_mean), width = 0.2) + 
  labs(y = 'Inward Bias (°)', x = 'Probe Distance (°)') +
  th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
             legend.position = 'none') +
  scale_y_continuous(breaks = seq(0, 25, 5)) + 
  scale_x_continuous(breaks = c(0, 30, 60, 90)) + 
  coord_capped_cart(ylim = c(-5, 25))
print(VS.grp)

ggsave(sprintf('VS_Group_%s.pdf', subDir), plot = VS.grp, height = 3, width = 4,  units = "in")

VS Quantile Analysis

VS.Quantile.plot <- VS.data %>%
  group_by(SN, Distance) %>%
  mutate(quantile = ntile(RT, 5)) %>%
  group_by(SN, Distance, quantile) %>%
  dplyr::summarise(bias_mean = mean(Bias, na.rm = TRUE), 
                   RT_mean = mean(RT, na.rm = TRUE)) %>%
  group_by(Distance, quantile) %>%
  dplyr::summarise(bias_grp_mean = mean(bias_mean, na.rm = TRUE), 
                   bias_std_err = std.error(bias_mean, na.rm = TRUE), 
                   RT_grp_mean = mean(RT_mean, na.rm = TRUE), 
                   RT_std_err = std.error(RT_mean, na.rm = TRUE)) %>%
  ggplot(aes(x = RT_grp_mean, y = bias_grp_mean, color = factor(Distance),group = factor(Distance))) + 
  geom_point() + 
  geom_line() + 
  geom_errorbarh(aes(xmax = RT_grp_mean + RT_std_err, xmin = RT_grp_mean - RT_std_err), alpha = 0.4)  +
  geom_errorbar(aes(ymax =bias_grp_mean + bias_std_err, ymin = bias_grp_mean - bias_std_err), alpha = 0.4) + 
  labs(y = 'Inward Bias (°)', x = 'RT (s)') +
  scale_color_manual(values = c("black", my.colors2[c(1, 2, 3)]), guide = FALSE) + 
  scale_x_continuous(breaks = seq(0.1, 0.5, 0.1)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
        legend.position = 'none') +
  th + 
  scale_y_continuous(breaks = c(0, 20, 40)) +
  coord_capped_cart(ylim = c(-5, 40),xlim = c(0.2, 0.5)) 

`summarise()` has grouped output by 'SN', 'Distance'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'Distance'. You can override using the `.groups` argument.

print(VS.Quantile.plot)

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 row(s) containing missing values (geom_path).
Warning: Removed 1 rows containing missing values (geom_errorbarh).
Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

ggsave(sprintf('VS_Quantile_%s.pdf', subDir), plot = VS.Quantile.plot, height = 3, width = 4,  units = "in")

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 row(s) containing missing values (geom_path).
Warning: Removed 1 rows containing missing values (geom_errorbarh).
Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

VS Quintile Normalize RT-Bias Dependency

VS.Quantile.norm.plot <- VS.data %>%
  group_by(SN, Distance) %>%
  filter(Distance != 0) %>%
  mutate(Scale_HandIB = scale(Hand_IB)) %>%
  mutate(quantile = ntile(RT, 5)) %>%
  group_by(SN, Distance, quantile) %>%
  dplyr::summarise(bias_mean = mean(Scale_HandIB, na.rm = TRUE), 
                   RT_mean = mean(RT, na.rm = TRUE)) %>%
  group_by(Distance, quantile) %>%
  dplyr::summarise(bias_grp_mean = mean(bias_mean, na.rm = TRUE), 
                   bias_std_err = std.error(bias_mean, na.rm = TRUE), 
                   RT_grp_mean = mean(RT_mean, na.rm = TRUE), 
                   RT_std_err = std.error(RT_mean, na.rm = TRUE)) %>%
  ggplot(aes(x = RT_grp_mean, y = bias_grp_mean, color = factor(Distance),group = factor(Distance))) + 
  geom_point() + 
  geom_line() + 
  geom_errorbarh(aes(xmax = RT_grp_mean + RT_std_err, xmin = RT_grp_mean - RT_std_err), alpha = 0.4)  +
  geom_errorbar(aes(ymax =bias_grp_mean + bias_std_err, ymin = bias_grp_mean - bias_std_err), alpha = 0.4) + 
  labs(y = 'Inward Bias (norm)', x = 'RT (s)') +
  scale_color_manual(values = c(my.colors2[c(1, 2, 3)]), guide = FALSE) + 
  scale_x_continuous(breaks = seq(0.1, 0.5, 0.1)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
        legend.position = 'none') +
  th + 
  coord_capped_cart(ylim = c(-1, 1),xlim = c(0.2, 0.5)) 

`summarise()` has grouped output by 'SN', 'Distance'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'Distance'. You can override using the `.groups` argument.

print(VS.Quantile.norm.plot)

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 row(s) containing missing values (geom_path).
Warning: Removed 1 rows containing missing values (geom_errorbarh).
Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

setwd("~/Desktop")

Warning: The working directory was changed to /Users/jonathantsay/Desktop inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

ggsave(sprintf('VS_Quantile_norm_%s.pdf', subDir), plot = VS.Quantile.norm.plot, height = 3, width = 4,  units = "in")

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 row(s) containing missing values (geom_path).
Warning: Removed 1 rows containing missing values (geom_errorbarh).
Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

Import Experiment 1

E1.data <- as_tibble( read.csv("/Users/jonathantsay/Dropbox/MOTOR/use_dependent/use-dependent-learning/UD_E1.csv", header = TRUE, sep = ',') %>%
                        dplyr::select(SN:RTb, hand_theta_40)) %>% mutate(group = 1) %>%
  mutate(Distance_raw = as.double( abs(ti - trainTgt)), 
         Distance = case_when(Distance_raw > 180 ~ abs(Distance_raw - 360), 
                              Distance_raw < -180 ~ abs(Distance_raw + 360),
                              Distance_raw <= 180 & Distance_raw >= -180 ~ Distance_raw),
         CN = TN, 
         Hand = hand_theta_40, Handb = NaN, RTb = NaN)

E1.data.base.out <- baseline_subtract(E1.data, 71, 140, 0)
E1.data <- E1.data.base.out[[1]]
E1.data.base <- E1.data.base.out[[2]]


E1.data <- E1.data %>%
  mutate(Hand_IB = case_when(trainTgt == 60 & (ti == 90   | ti == 120  | ti == 150 ) ~ -Handb, 
                             trainTgt == 60 & (ti == 0   | ti == 30  | ti == 330 | ti == 60) ~ Handb, 
                             trainTgt == 150 & (ti == 180  | ti == 210  | ti == 240) ~ -Handb, 
                             trainTgt == 150 & (ti == 60   | ti == 90  | ti == 120 | ti == 150) ~ Handb))


  

Experiment 1 Individual Data

for(si in unique(E1.data$SN)){
  
  groupvar <- unique( E1.data$group[E1.data$SN == si] ) 
  myymin <- min( E1.data$RT[E1.data$SN == si & E1.data$block >= 3 & E1.data$fbi == 0] ) + 0.1
  myymax <- max( E1.data$RT[E1.data$SN == si & E1.data$block >= 3 & E1.data$fbi == 0] ) - 0.1
  
  ind.plot <- E1.data %>% 
    filter(SN == si & Distance != 0) %>%  
    filter(block >= 3 & fbi == 0 ) %>%
    ggplot(aes(x = RT, y = Hand_IB, group = factor(Distance))) + 
    geom_point(alpha = 0.4, aes(color = factor(Distance))) + 
    geom_smooth(method = 'lm', aes(color = factor(Distance)), size = 0.5, se = FALSE, linetype= groupvar) +
    #stat_cor(method = 'pearson')+ 
    facet_rep_wrap(.~Distance, repeat.tick.labels = TRUE, ncol = 4) + 
    scale_x_continuous(breaks = c(round(myymin, 1), round(myymax, 1) )) +
    scale_color_manual(values = c(my.colors2[c(1, 2, 3)]), guide = FALSE) + 
    labs(y = 'Inward Bias (°)', x = 'RT (s)', subtitle = sprintf('Exp 1 Subject #%s', si)) +
    th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
               legend.position = 'none') +
    scale_y_continuous(breaks = seq(-50, 100, 50)) + 
    coord_capped_cart(ylim = c(-50, 100)) 
  print(ind.plot) 
  
  ggsave(sprintf('E1_SUB%s_Grp%s_notitle_%s.pdf', si, groupvar, subDir), plot = ind.plot, height = 3, width = 4,  units = "in")
}

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

Experiment 1 Group Data

E1.ind <- E1.data %>%
  filter(block >= 3 & fbi == 0 ) %>%
  group_by(SN, Distance, group) %>%
  dplyr::summarise(hand_mean = mean(Hand_IB, na.rm = TRUE), 
                   rt_mean = mean(RT, na.rm = TRUE))

`summarise()` has grouped output by 'SN', 'Distance'. You can override using the `.groups` argument.

E1.ind.base <- E1.data %>%
  filter(block< 3 & fbi == 0 ) %>%
  group_by(SN, Distance, group) %>%
  dplyr::summarise(hand_mean = mean(Hand_IB, na.rm = TRUE), 
                   rt_mean = mean(RT, na.rm = TRUE))

`summarise()` has grouped output by 'SN', 'Distance'. You can override using the `.groups` argument.

E1.grp <- E1.ind %>% 
  group_by(Distance, group) %>%
  dplyr::summarise(hand_grp_mean = mean(hand_mean, na.rm = TRUE),
                   hand_grp_stderr = std.error(hand_mean, na.rm = TRUE),
                   rt_grp_mean = mean(rt_mean, na.rm = TRUE), 
                   rt_grp_stderr = std.error(rt_mean, na.rm = TRUE))

`summarise()` has grouped output by 'Distance'. You can override using the `.groups` argument.

E1.grp.plot <- E1.ind %>%
  ggplot(aes(x = Distance, y = hand_mean, group = group)) + 
  geom_segment(aes(x=-5, xend=95, y = 0, yend = 0), color = 'lightgrey') +
  stat_summary(fun = 'mean', geom = "point") + 
  stat_summary(fun = 'mean', geom = "line", aes(linetype = factor(group))) + 
  stat_summary(fun.data = 'mean_se', geom = "errorbar", width = 2) + 
  labs(y = 'Inward Bias (°)', x = 'Probe Distance (°)') +
  scale_x_continuous(breaks = c(0, 30, 60, 90)) +
  th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
             legend.position = 'none') +
  coord_capped_cart(ylim = c(-5, 40))
print(E1.grp.plot)  

ggsave(sprintf('E1_Group_%s.pdf', subDir), plot = E1.grp.plot, height = 3, width = 4,  units = "in")

Experiment 1 Group Data, with First Peak (Obtained from MClust)

firstpeak <- setNames ( as.data.frame(matrix(nrow = 4, ncol = 5)), 
                        c("Distance", "group", "hand_mean", "hand_ub", "hand_lb") )

firstpeak$Distance <- c(0, 30, 60, 90)
firstpeak$hand_mean <- c(-1.1, 1.1, 4.2, 3.4)
firstpeak$hand_ub <- c(-1.1, 2.3, 5.1, 4.2)
firstpeak$hand_lb <- c(-1.1, 0.1, 3.2, 2.5)

E1.grpNoJit.plot <- E1.ind %>%
  filter(group == 1) %>%
  ggplot(aes(x = Distance, y = hand_mean, group = group)) + 
  geom_segment(aes(x=-5, xend=95, y = 0, yend = 0), color = 'lightgrey') +
  stat_summary(fun = 'mean', geom = "point") + 
  stat_summary(fun = 'mean', geom = "line", aes(linetype = factor(group))) + 
  stat_summary(fun.data = 'mean_se', geom = "errorbar", width = 2) + 
  geom_point(inherit.aes = FALSE, data = firstpeak, aes(x = Distance, y = hand_mean), color = "black")+
  geom_errorbar(inherit.aes = FALSE, data = firstpeak, aes(x = Distance, ymin = hand_lb, ymax = hand_ub), color = "black", width = 2)+
  geom_line(inherit.aes = FALSE, data = firstpeak, aes(x = Distance, y = hand_mean), color = "black", linetype = "dashed")+
  labs(y = 'Inward Bias (°)', x = 'Probe Distance (°)') +
  scale_x_continuous(breaks = c(0, 30, 60, 90)) +
  th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
             legend.position = 'none') +
  coord_capped_cart(ylim = c(-5, 40))
print(E1.grpNoJit.plot)  

ggsave(sprintf('E1_GroupNoJit_%s.pdf', subDir), plot = E1.grpNoJit.plot, height = 3, width = 4,  units = "in")

std.error(E1.ind$hand_mean[E1.ind$group == 1 & E1.ind$Distance == 0])

[1] 1.012352

Experiment 1 Quantile Analysis

E1.ind.Quantile.data <- E1.data %>%
  filter(block >= 3 & fbi == 0 ) %>%
  group_by(SN, Distance, group) %>%
  mutate(quantile = ntile(RT, 5)) %>%
  group_by(SN, Distance, quantile, group) %>%
  dplyr::summarise(bias_mean = mean(Hand_IB, na.rm = TRUE), 
                   RT_mean = mean(RT, na.rm = TRUE))

`summarise()` has grouped output by 'SN', 'Distance', 'quantile'. You can override using the `.groups` argument.

E1.Quantile.data <- E1.ind.Quantile.data %>% 
  group_by(Distance, quantile, group) %>%
  dplyr::summarise(bias_grp_mean = mean(bias_mean, na.rm = TRUE), 
                   bias_std_err = std.error(bias_mean, na.rm = TRUE), 
                   RT_grp_mean = mean(RT_mean, na.rm = TRUE), 
                   RT_std_err = std.error(RT_mean, na.rm = TRUE))

`summarise()` has grouped output by 'Distance', 'quantile'. You can override using the `.groups` argument.

E1.Quantile.grp1.plot <- E1.Quantile.data %>%
  filter(group == 1) %>%
  ggplot(aes(x = RT_grp_mean, y = bias_grp_mean, color = factor(Distance), group = factor(Distance))) + 
  geom_point() + 
  geom_line() + 
  geom_errorbarh(aes(xmax = RT_grp_mean + RT_std_err, xmin = RT_grp_mean - RT_std_err), alpha = 0.4)  +
  geom_errorbar(aes(ymax =bias_grp_mean + bias_std_err, ymin = bias_grp_mean - bias_std_err), alpha = 0.4) + 
  labs(y = 'Inward Bias (°)', x = 'RT (s)') +
  scale_color_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) + 
  scale_x_continuous(breaks = seq(0.2, 0.6, 0.1)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
        legend.position = 'none') +
  th + 
  scale_y_continuous(breaks = c(0, 25, 50, 75)) +
  coord_capped_cart(ylim = c(-5, 75), xlim = c(0.175, 0.6)) 
print(E1.Quantile.grp1.plot)

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

ggsave(sprintf('E1_Quantile_Grp1_%s.pdf', subDir), plot = E1.Quantile.grp1.plot, height = 3, width = 4,  units = "in")

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

Experiment 1 Normalized Quintile Analysis

E1.ind.Quantile.norm.data <- E1.data %>%
  filter(block >= 3 & fbi == 0 & group == 1) %>%
  group_by(SN, Distance) %>%
  mutate(Scale_HandIB = scale(Hand_IB)) %>%
  mutate(quantile = ntile(RT, 5)) %>%
  group_by(SN, Distance, quantile) %>%
  dplyr::summarise(bias_mean = mean(Scale_HandIB, na.rm = TRUE), 
                   RT_mean = mean(RT, na.rm = TRUE))

`summarise()` has grouped output by 'SN', 'Distance'. You can override using the `.groups` argument.

E1.Quantile.norm.data <- E1.ind.Quantile.norm.data %>% 
  group_by(Distance, quantile) %>%
  dplyr::summarise(bias_grp_mean = mean(bias_mean, na.rm = TRUE), 
                   bias_std_err = std.error(bias_mean, na.rm = TRUE), 
                   RT_grp_mean = mean(RT_mean, na.rm = TRUE), 
                   RT_std_err = std.error(RT_mean, na.rm = TRUE))

`summarise()` has grouped output by 'Distance'. You can override using the `.groups` argument.

E1.Quantile.all.norm.plot <- E1.Quantile.norm.data %>%
  filter(Distance != 0) %>%
  ggplot(aes(x = RT_grp_mean, y = bias_grp_mean, color = factor(Distance), group = factor(Distance))) + 
  geom_point() + 
  geom_line() + 
  geom_errorbarh(aes(xmax = RT_grp_mean + RT_std_err, xmin = RT_grp_mean - RT_std_err), alpha = 0.4)  +
  geom_errorbar(aes(ymax =bias_grp_mean + bias_std_err, ymin = bias_grp_mean - bias_std_err), alpha = 0.4) + 
  labs(y = 'Inward Bias (norm)', x = 'RT (s)') +
  scale_color_manual(values = c(my.colors2[c(1, 2, 3)]), guide = FALSE) + 
  scale_x_continuous(breaks = seq(0.2, 0.7, 0.1)) + 
    th + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
        legend.position = 'none') 
print(E1.Quantile.all.norm.plot)

ggsave(sprintf('E1_Quantile_Grp1_Norm_%s.pdf', subDir), plot = E1.Quantile.all.norm.plot, height = 3, width = 4,  units = "in")

Trajectory

E1.traj.data <- as_tibble( read.csv("/Users/jonathantsay/Dropbox/MOTOR/use_dependent/UD_WithTraj.csv", header = TRUE, sep = ',')) %>%
  dplyr::select(SN:hx_138, hy_1:hy_138) 


for(si in unique(E1.traj.data$SN)){
  E1.traj.data.long <- E1.traj.data %>% 
    filter(SN == si) %>%
    gather(Cond, Value, hx_1:hy_138) %>%
    separate(Cond, into = c('XorY', 'Samplepoint', sep = "_")) %>%
    spread(XorY, Value)
  
  ind.traj.plot <- E1.traj.data.long %>%
    filter(block >= 3, fbi == 0) %>%
    drop_na(hx, hy) %>%  
    ggplot(aes(x = hx, y =hy, color = factor(ti), group = TN)) + 
    geom_path(size = 0.5) +
    coord_capped_cart(ylim = c(-150, 150), xlim = c(-150, 150)) +
    th +
    theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
          legend.position = 'none') +
    labs(x = 'x (mm)', y = 'y (mm)')
  
  print(ind.traj.plot)
  
  ggsave(sprintf('E1_SUB%s_Traj_%s.pdf', si, subDir), plot = ind.traj.plot, height = 3, width = 4,  units = "in")
  
}

Load in E2 (Reward exp) Hold time 500, with forced delay.

E2.data <- as_tibble( read.csv("/Users/jonathantsay/Dropbox/MOTOR/use_dependent/use-dependent-learning/UD_E2.csv", header = TRUE, sep = ',')) %>% 
  mutate(Distance_raw = as.double( abs(ti - trainTgt)), Distance = case_when(Distance_raw > 180 ~ abs(Distance_raw - 360),
                                                    Distance_raw < -180 ~ abs(Distance_raw + 360),
                                                    Distance_raw <= 180 & Distance_raw >= -180 ~ Distance_raw),
                                                    CN = TN,
                                                    Hand = hand_theta_40, Handb = NaN, RTb = NaN)


E2.data.base.out <- baseline_subtract(E2.data, 71, 140, 0)
E2.data <- E2.data.base.out[[1]]
E2.data.base <- E2.data.base.out[[2]]

# flip to inward bias 
E2.data <- E2.data %>%
  mutate(Hand_IB = case_when(trainTgt == 60 & (ti == 90   | ti == 120  | ti == 150 ) ~ -Handb, 
                             trainTgt == 60 & (ti == 0   | ti == 30  | ti == 330 | ti == 60) ~ Handb, 
                             trainTgt == 150 & (ti == 180  | ti == 210  | ti == 240) ~ -Handb, 
                             trainTgt == 150 & (ti == 60   | ti == 90  | ti == 120 | ti == 150) ~ Handb))

Experiment 2 Group Data

E2.ind <- E2.data %>%
  filter(block >= 3 & fbi == 0 ) %>%
  group_by(SN, Distance, group) %>%
  dplyr::summarise(hand_mean = mean(Hand_IB, na.rm = TRUE), 
                   rt_mean = mean(RT, na.rm = TRUE), 
                   mt_mean = mean(MT, na.rm = TRUE))

`summarise()` has grouped output by 'SN', 'Distance'. You can override using the `.groups` argument.

E2.grp <- E2.ind %>% 
  group_by(Distance, group) %>%
  dplyr::summarise(hand_grp_mean = mean(hand_mean, na.rm = TRUE),
                   hand_grp_stderr = std.error(hand_mean, na.rm = TRUE),
                   rt_grp_mean = mean(rt_mean, na.rm = TRUE), 
                   rt_grp_stderr = std.error(rt_mean, na.rm = TRUE))

`summarise()` has grouped output by 'Distance'. You can override using the `.groups` argument.

E2.grpBoth.plot <- E2.ind %>%
  ggplot(aes(x = Distance, y = hand_mean, group = group, color = group)) + 
  geom_segment(aes(x=-5, xend=95, y = 0, yend = 0), color = 'lightgrey') +
  stat_summary(fun = 'mean', geom = "point") + 
  stat_summary(fun = 'mean', geom = "line", linetype = 1) + 
  stat_summary(fun.data = 'mean_se', geom = "errorbar", width = 2) + 
  labs(y = 'Inward Bias (°)', x = 'Probe Distance (°)') +
  scale_x_continuous(breaks = c(0, 30, 60, 90)) +
  scale_color_manual(values = c( 'darkgrey', 'black')) +
  th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
             legend.position = 'none') +
  coord_capped_cart(ylim = c(-5, 10))
print(E2.grpBoth.plot)  

ggsave(sprintf('E2_GroupBoth_%s.pdf', subDir), plot = E2.grpBoth.plot, height = 3, width = 4,  units = "in")

Experiment 2: Quantile Analysis Combined

E2.data <- E2.data %>% mutate(SN_new = ifelse(group == "R", SN + 999, SN))

E2.Quantile.comb.data <-E2.data %>%
  filter(block >= 3 & fbi == 0 ) %>%
  group_by(SN_new, Distance) %>%
  mutate(quantile = ntile(RT, 5)) %>%
  group_by(SN_new, Distance, quantile) %>%
  dplyr::summarise(bias_mean = mean(Hand_IB, na.rm = TRUE), 
                   RT_mean = mean(RT, na.rm = TRUE)) %>%
  group_by(Distance, quantile) %>%
  dplyr::summarise(bias_grp_mean = mean(bias_mean, na.rm = TRUE), 
                   bias_std_err = std.error(bias_mean, na.rm = TRUE), 
                   RT_grp_mean = mean(RT_mean, na.rm = TRUE), 
                   RT_std_err = std.error(RT_mean, na.rm = TRUE))

`summarise()` has grouped output by 'SN_new', 'Distance'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'Distance'. You can override using the `.groups` argument.

E2.Quantile.comb.plot <- E2.Quantile.comb.data %>%
  ggplot(aes(x = RT_grp_mean, y = bias_grp_mean, color = factor(Distance), group = factor(Distance))) + 
  geom_point() + 
  geom_line() + 
  geom_errorbarh(aes(xmax = RT_grp_mean + RT_std_err, xmin = RT_grp_mean - RT_std_err), alpha = 0.4)  +
  geom_errorbar(aes(ymax =bias_grp_mean + bias_std_err, ymin = bias_grp_mean - bias_std_err), alpha = 0.4) + 
  labs(y = 'Inward Bias (°)', x = 'RT (s)') +
  scale_color_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) + 
  #scale_x_continuous(breaks = seq(0.2, 0.7, 0.1)) + 
  scale_y_continuous(breaks = seq(0, 10, 5)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
        legend.position = 'none') +
  th + 
  coord_capped_cart(ylim = c(-2, 10)) 
print(E2.Quantile.comb.plot)

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

ggsave(sprintf('E2_Quantile_Comb_%s.pdf', subDir), plot = E2.Quantile.comb.plot, height = 3, width = 4,  units = "in")

Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

for(si in unique(E2.data$SN)){
  
  pro.data <- E2.data %>% 
    filter(SN == si) %>%  
    filter(block >= 3 & fbi == 0 & Distance > 0 & RT < 3)
  
  groupvar <- 1
  myymin <- min( pro.data$RT[pro.data$SN == si & pro.data$block >= 3 & pro.data$fbi == 0], na.rm = TRUE) + 0.2
  myymax <- max( pro.data$RT[pro.data$SN == si & pro.data$block >= 3 & pro.data$fbi == 0], na.rm = TRUE) - 0.2
  
  ind.plot <- pro.data %>%
    ggplot(aes(x = RT, y = Hand_IB, group = factor(Distance))) + 
    geom_point(alpha = 0.4, aes(color = factor(Distance))) + 
    geom_smooth(method = 'lm', aes(color = factor(Distance)), size = 0.5, se = FALSE, linetype= groupvar) + 
    facet_rep_wrap(.~Distance, repeat.tick.labels = TRUE, ncol = 4) + 
    #scale_x_continuous(breaks = c(round(myymin, 1), round(myymax, 1) )) +
    scale_color_manual(values = c(my.colors2[c(1, 2, 3)]), guide = FALSE) + 
    labs(y = 'Inward Bias (°)', x = 'RT (s)', subtitle = sprintf('Exp 2 Subject #%s', si)) + 
    th + theme(axis.text.x = element_text(angle = 45, hjust = 0.95), 
               legend.position = 'none') +
    coord_capped_cart(ylim = c(-25, 25)) 
  print(ind.plot) 
  
  ggsave(sprintf('aE2_SUB%s_Grp%s_%s.pdf', si, groupvar, subDir), plot = ind.plot, height = 3, width = 4,  units = "in")
  
}

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

rewardcdf <- E2.data %>%
  filter(block >= 3 & Distance == 0) %>%
  ggplot(aes(x = abs(hand_theta), alpha = group)) + 
  stat_ecdf(geom = "step") + 
  scale_alpha_manual(values = c(0.2, 1), guide = FALSE) +
  scale_x_continuous(limits = c(0, 40)) +
  labs(y = 'Cumulative Probability', x = 'Error (°)') +
  th
print(rewardcdf)

Warning: Removed 13 rows containing non-finite values (stat_ecdf).
Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

ggsave(sprintf('E2_rewardcdf_%s.pdf', subDir), plot = rewardcdf, height = 3, width = 3,  units = "in")

Warning: Removed 13 rows containing non-finite values (stat_ecdf).
Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

E1 Bimodal

E1.data$Distance <- factor(E1.data$Distance, levels = c("0", "30", "60", "90"))
E1.data.density <- E1.data %>%
  filter(block > 2) %>%
  ggplot(aes(x = Hand_IB, fill = Distance)) +
  geom_density(alpha = 0.5) +
  scale_x_continuous(breaks = seq(-30, 120, 30)) +
  facet_rep_wrap(Distance~., ncol = 4, repeat.tick.labels = TRUE) + 
  th + theme(legend.position = "none") + 
  scale_fill_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) +
  coord_capped_cart(xlim = c(-30, 120)) +
  labs(x = "Heading Angle (°)", y = "Probability Density")
print(E1.data.density)

#setwd("~/Desktop")
ggsave(sprintf('E1_Density_%s.pdf', subDir), plot = E1.data.density, height = 3, width = 10.5,  units = "in")

VS.data.allprobes <-  as_tibble( read.csv("/Users/jonathantsay/Dropbox/MOTOR/use_dependent/use-dependent-learning/VerstynenSabes2011.csv", header = TRUE, sep = ',')) %>% mutate(Bias = Hand_IB)

VS.data.density <- VS.data.allprobes %>%
  ggplot(aes(x = Hand_IB, fill = factor(Distance)))  +
  geom_density(alpha = 0.5) +
   scale_x_continuous(breaks = seq(-30, 120, 30), limits = c(-60, 150)) +
  facet_rep_wrap(Distance~., ncol = 4, repeat.tick.labels = TRUE) + 
  th + theme(legend.position = "none") + 
  scale_fill_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) +
  labs(x = "Heading Angle (°)", y = "Probability Density")
print(VS.data.density)

Warning: Removed 1 rows containing non-finite values (stat_density).

setwd("~/Desktop")

Warning: The working directory was changed to /Users/jonathantsay/Desktop inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

#ggsave(sprintf('VS_Density_%s.pdf', subDir), plot = VS.data.density, height = 3, width = 10.5,  units = "in")

E2 Bimodal

E2.data$Distance <- factor(E2.data$Distance, levels = c("0", "30", "60", "90"))
E2.data.density <- E2.data %>%
  filter(block > 2) %>%
  ggplot(aes(x = Hand_IB, fill = Distance)) +
  #geom_histogram(binwidth = 5) +
  geom_density(alpha = 0.5) +
  #geom_vline(xintercept = 90, linetype = "dotted") +
  scale_x_continuous(breaks = seq(-30, 120, 30), limits = c(-60, 150)) +
  coord_capped_cart(xlim = c(-30,120)) +
  facet_rep_wrap(Distance~., ncol = 4, repeat.tick.labels = TRUE) + 
  th + theme(legend.position = "none") + 
  scale_fill_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) +
  labs(x = "Heading Angle (°)", y = "Probability Density")
print(E2.data.density)

Warning: Removed 10 rows containing non-finite values (stat_density).

#setwd("~/Desktop")
ggsave(sprintf('E2_Density_%s.pdf', subDir), plot = E2.data.density, height = 3, width = 10.5,  units = "in")

Warning: Removed 10 rows containing non-finite values (stat_density).

E1 modeling

E1.data$Distance <- factor(E1.data$Distance, levels = c(0, 90, 30, 60))
E1.data.mod <- E1.data %>% filter(group == 1 & block >= 3) %>% group_by(SN, Distance) %>% mutate(scale_Hand_IB = scale(Hand_IB))
myE1mod <- lmer(scale_Hand_IB ~ RT * factor(Distance) + (1 |SN), data = E1.data.mod %>% filter(Distance != "0"))

boundary (singular) fit: see ?isSingular
Warning: Model failed to converge with 1 negative eigenvalue: -1.5e+04

summary(myE1mod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale_Hand_IB ~ RT * factor(Distance) + (1 | SN)
   Data: E1.data.mod %>% filter(Distance != "0")

REML criterion at convergence: 2581.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7395 -0.6803 -0.0354  0.6034  3.7991 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 SN       (Intercept) 4.527e-34 2.128e-17
 Residual             8.550e-01 9.247e-01
Number of obs: 960, groups:  SN, 10

Fixed effects:
                        Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)             0.675877   0.115995 954.000000   5.827 7.72e-09 ***
RT                     -1.808606   0.277871 954.000000  -6.509 1.22e-10 ***
factor(Distance)30      0.177712   0.170048 954.000000   1.045   0.2963    
factor(Distance)60     -0.007582   0.170643 954.000000  -0.044   0.9646    
RT:factor(Distance)30  -0.854560   0.449122 954.000000  -1.903   0.0574 .  
RT:factor(Distance)60  -0.123912   0.431110 954.000000  -0.287   0.7738    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) RT     f(D)30 f(D)60 RT:(D)3
RT          -0.895                             
fctr(Dst)30 -0.682  0.611                      
fctr(Dst)60 -0.680  0.609  0.464               
RT:fct(D)30  0.554 -0.619 -0.900 -0.376        
RT:fct(D)60  0.577 -0.645 -0.394 -0.903  0.399 
convergence code: 0
boundary (singular) fit: see ?isSingular

Anova(myE1mod, type = 'III')

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: scale_Hand_IB
                      Chisq Df Pr(>Chisq)    
(Intercept)         33.9515  1  5.650e-09 ***
RT                  42.3644  1  7.576e-11 ***
factor(Distance)     1.4487  2     0.4846    
RT:factor(Distance)  3.8846  2     0.1434    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

eta_sq(myE1mod)

Warning: Model failed to converge with 1 negative eigenvalue: -1.5e+04
Warning: Currently only supports partial eta squared for this class of objects.

                 term etasq
1                  RT 0.122
2    factor(Distance) 0.002
3 RT:factor(Distance) 0.004

VS modeling

VS.data$Distance <- factor(VS.data$Distance, levels = c(0, 30, 60, 90))
VS.data <- VS.data %>% group_by(SN, Distance) %>% mutate(scale_hand_IB = scale(Hand_IB))
myVSmod <- lmer(Hand_IB ~  RT * factor(Distance) + (1|SN), data = VS.data)
summary(myVSmod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hand_IB ~ RT * factor(Distance) + (1 | SN)
   Data: VS.data

REML criterion at convergence: 5506.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6100 -0.5525 -0.0734  0.4984  4.1101 

Random effects:
 Groups   Name        Variance Std.Dev.
 SN       (Intercept)  43.62    6.604  
 Residual             223.44   14.948  
Number of obs: 671, groups:  SN, 8

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              4.054      7.746  335.767   0.523  0.60108    
RT                     -12.227     22.841  660.334  -0.535  0.59261    
factor(Distance)30      29.856      9.208  657.002   3.242  0.00125 ** 
factor(Distance)60      57.332      8.866  656.755   6.466 1.96e-10 ***
factor(Distance)90      97.954      8.775  656.919  11.163  < 2e-16 ***
RT:factor(Distance)30  -67.471     28.245  656.910  -2.389  0.01719 *  
RT:factor(Distance)60 -137.197     26.567  656.749  -5.164 3.21e-07 ***
RT:factor(Distance)90 -232.285     26.214  656.894  -8.861  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) RT     f(D)30 f(D)60 f(D)90 RT:(D)3 RT:(D)6
RT          -0.933                                            
fctr(Dst)30 -0.639  0.650                                     
fctr(Dst)60 -0.687  0.701  0.581                              
fctr(Dst)90 -0.694  0.707  0.586  0.615                       
RT:fct(D)30  0.627 -0.672 -0.979 -0.569 -0.575                
RT:fct(D)60  0.698 -0.748 -0.581 -0.976 -0.617  0.598         
RT:fct(D)90  0.708 -0.759 -0.588 -0.618 -0.976  0.606   0.652 

Anova(myVSmod, type = 'III')

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: Hand_IB
                       Chisq Df Pr(>Chisq)    
(Intercept)           0.2739  1     0.6007    
RT                    0.2866  1     0.5924    
factor(Distance)    142.1544  3     <2e-16 ***
RT:factor(Distance)  92.7832  3     <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

eta_sq(myVSmod)

Warning: Currently only supports partial eta squared for this class of objects.

                 term etasq
1                  RT 0.150
2    factor(Distance) 0.178
3 RT:factor(Distance) 0.124

E2 modeling

E2.data$Distance <- factor(E2.data$Distance, levels = c(0, 30, 60, 90))
myE2mod <- lmer(Hand_IB ~ RT*Distance + (SN_new| group), data = E2.data %>% filter(block >= 3 & fbi == 0 ))

boundary (singular) fit: see ?isSingular

summary(myE2mod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Hand_IB ~ RT * Distance + (SN_new | group)
   Data: E2.data %>% filter(block >= 3 & fbi == 0)

REML criterion at convergence: 178363.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-20.8926  -0.5479   0.0148   0.5748  19.9080 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr 
 group    (Intercept) 1.316e+02 11.4730       
          SN_new      8.761e-03  0.0936  -1.00
 Residual             8.304e+01  9.1129       
Number of obs: 24576, groups:  group, 2

Fixed effects:
                Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)      24.8318     1.5932     7.3262  15.586 6.96e-07 ***
RT                3.8310     0.2148 24523.8910  17.836  < 2e-16 ***
Distance30        2.4831     1.2907 24560.4302   1.924   0.0544 .  
Distance60        6.3609     1.5142 24560.9518   4.201 2.67e-05 ***
Distance90        7.7762     1.0511 24559.8088   7.398 1.42e-13 ***
RT:Distance30     0.6921     1.3796 24560.4580   0.502   0.6159    
RT:Distance60    -2.8172     1.6541 24561.0021  -1.703   0.0886 .  
RT:Distance90    -5.4704     1.1124 24559.7931  -4.918 8.82e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) RT     Dstn30 Dstn60 Dstn90 RT:D30 RT:D60
RT          -0.183                                          
Distance30  -0.007  0.152                                   
Distance60  -0.002  0.129  0.022                            
Distance90  -0.020  0.187  0.032  0.027                     
RT:Distnc30  0.005 -0.154 -0.977 -0.021 -0.029              
RT:Distnc60  0.000 -0.128 -0.020 -0.983 -0.024  0.020       
RT:Distnc90  0.018 -0.192 -0.030 -0.025 -0.965  0.030  0.025
convergence code: 0
boundary (singular) fit: see ?isSingular

Anova(myE2mod, type = 'III')

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: Hand_IB
              Chisq Df Pr(>Chisq)    
(Intercept) 242.913  1  < 2.2e-16 ***
RT          318.134  1  < 2.2e-16 ***
Distance     73.326  3  8.275e-16 ***
RT:Distance  27.143  3  5.494e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

eta_sq(myE2mod)

Warning: Currently only supports partial eta squared for this class of objects.

         term etasq
1          RT 0.000
2    Distance 0.003
3 RT:Distance 0.001

Make overlapping distribution

E1.data$Distance <- factor(E1.data$Distance, levels = c("0", "30", "60", "90"))
E1.data.density.sub <- E1.data %>%
  filter(block > 2) %>%
  ggplot(aes(x = Hand_IB, fill = Distance, group = factor(SN))) +
  #geom_histogram(binwidth = 5) +
  #geom_vline(xintercept = 90, linetype = "dotted") +
  geom_density(alpha = 0.1, size = 0.1) +
  scale_x_continuous(breaks = seq(-30, 120, 30)) +
  facet_rep_wrap(Distance~., ncol = 4, repeat.tick.labels = TRUE) + 
  th + theme(legend.position = "none") + 
  scale_fill_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) +
  coord_capped_cart(xlim = c(-30, 120)) +
  labs(x = "Heading Angle (°)", y = "Probability Density")
print(E1.data.density.sub)

ggsave(sprintf('E1_ind_density_%s.pdf', subDir), plot = E1.data.density.sub, height = 3, width = 10.5,  units = "in")

E2.data$Distance <- factor(E2.data$Distance, levels = c("0", "30", "60", "90"))
E2.data.density.sub <- E2.data %>%
  filter(block > 2) %>%
  ggplot(aes(x = Hand_IB, fill = Distance, group = factor(SN))) +
  #geom_histogram(binwidth = 5) +
  geom_density(alpha = 0.1, size = 0.1) +
  #geom_vline(xintercept = 90, linetype = "dotted") +
  scale_x_continuous(breaks = seq(-30, 120, 30), limits = c(-60, 150)) +
  coord_capped_cart(xlim = c(-30,120)) +
  facet_rep_wrap(Distance~., ncol = 4, repeat.tick.labels = TRUE) + 
  th + theme(legend.position = "none") + 
  scale_fill_manual(values = c('black', my.colors2[c(1, 2, 3)]), guide = FALSE) +
  labs(x = "Heading Angle (°)", y = "Probability Density")
print(E2.data.density.sub)

Warning: Removed 10 rows containing non-finite values (stat_density).

ggsave(sprintf('E2_ind_density_%s.pdf', subDir), plot = E2.data.density.sub, height = 3, width = 10.5,  units = "in")

Warning: Removed 10 rows containing non-finite values (stat_density).